Equilibrium-like behavior in far-from-equilibrium chemical reaction networks 
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In an equilibrium chemical reaction mixture, the number of molecules present obeys a Poisson 
distribution. We ask when the same is true of the steady state of a nonequilibrium reaction network 
and obtain an essentially complete answer. In particular, we show that networks with certain topo- 
logical features must have a Poisson distribution, whatever the reaction rates. Such driven systems 
also obey an analog of the fluctuation-dissipation theorem. Our results may be relevant to biologi- 
cal systems and to the larger question of how equilibrium concepts might apply to nonequilibrium 
systems. 



Physicists have long sought to extend the framework 
of equilibrium statistical mechanics to describe systems 
far from equilibrium. Although this general goal is still 
far from attained, recent years have seen the discovery 
of several specific systems that, while strongly driven, 
nonetheless act very much as if they were in equilib- 
rium. Such unexpected thermal behavior has been ob- 
served, for example, in groups of fluidized spheres in 
vibrated granular matter [2|, and in slowly sheared col- 
lections of particles near the jamming transition [H, 0]. 
These cases are striking for the extent to which an equi- 
librium analogy holds; the dynamics of fluctuations and 
sometimes entire distributions of observables correspond 
to those expected in equilibrium. Such observations serve 
as tantalizing hints that a simpler structure may be hid- 
ing under the apparent complexities of nonequilibrium 
statistical physics. It is thus of considerable interest to 
understand why these examples act as they do and what 
the limits are of their likeness to thermal systems. In 
the cases just cited, the equilibrium-like properties were 
found either experimentally or through large-scale simu- 
lations, and a full theoretical understanding of the obser- 
vations has been slow to develop. In particular, despite a 
few important contributions [5] , it has proven difficult to 
say precisely when an equilibrium-like description should 
hold. Here, we address this question for a new class of 
systems, chemical reaction networks, and show that in 
this case it can be answered in some detail. 

The statistical physics of chemical systems has re- 
cently benefited from a resurgence of interest motivated 
by biological applications Q. In some cells, regulatory 
molecules can be present in as few as one or two copies, 
on average. In this regime, the intrinsic stochasticity of 
chemical reactions causes large fluctuations in molecule 
numbers, and it is natural to ask how cells can continue 
to function in the face of the resulting uncertainty. A 
complete answer to this question requires knowing what 
determines the molecule number distribution in driven, 
reacting systems. We take an important step in this di- 
rection by giving a full characterization of those reac- 
tion networks whose steady-state distribution (SSD) of 
molecule numbers is Poissonian, as it would be in an 
equilibrium system. We show that a surprisingly large 
class of networks have this property. Such networks can 
include completely irreversible reactions and an arbitrar- 



ily large number of chemical species; the network topol- 
ogy should, however, satisfy a sparseness condition that 
allows the construction of a generalized free energy land- 
scape. 

In what follows, we first show that whether a mass- 
action network has a Poisson SSD (PSSD) is closely re- 
lated to whether the mean-field kinetic equations for the 
same network have a unique fixed point. This latter 
problem has a long history in the chemical engineering 
literature @, [1] , and by adapting some results from the 
deterministic case, we are able to give a mathematical 
description of the set of networks with a PSSD. This 
set includes, for any choice of reaction rates, networks 
with a so-called deficiency zero topology. We then show 
that the equilibrium analogy extends beyond the form 
of the SSD. Networks with a PSSD satisfy a version of 
the fluctuation-dissipation theorem, and when two such 
systems are brought into contact, they exchange parti- 
cles to equalize their values of a generalized chemical 
potential. Analogous results hold for certain non-mass- 
action kinetic schemes. A significant group of far-from- 
equilibrium networks thus shows a marked resemblance 
to equilibrium systems, and it is possible to gain consid- 
erable insight into when and why this occurs. 

We begin with some formal preliminaries [7|, [8[. By 
a chemical reaction network, we mean a set of chemical 
species (or molecules) that can undergo a prescribed set 
of reactions (Fig.[TJ. For most of this paper, we consider 
a well-stirred reaction vessel, so that spatial degrees of 
freedom are irrelevant; later, we will argue that our re- 
sults extend to models that explicitly include molecular 
diffusion. We denote species by uppercase letters. A 
complex is the collection of species on one side of a reac- 
tion arrow together with their stoichiometric coefficients. 
In the reaction 2 A + B — ► C, the ng — 3 species are 
{A, B, C}, and the n c = 2 complexes are {2A + B, C}. A 
network may include pseudo-reactions that, for example, 
absorb the concentration of a reactant that is in excess 
into the rate constant; thus, a reaction like <f> — > A, in- 
dicating that A molecules are created out of nothing, is 
possible. 

Chemical reactions in dilute systems are most com- 
monly modeled by deterministic rate equations (DRE). 
In this description, the relevant variables are the 
concentrations c = (cai, ca 2) ca 3 > • • •) of the species 
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FIG. 1: Some chemical reaction networks. a)-c) have 
complexes {A, 2B, C, 3D} and I = 1 linkage classes. d) 
has complexes {A, 2B, A+C, D, B+E} and 1 = 2. e) 
has complexes {A+E, C, B+E} and I = 1, f) complexes 
{A+E, C, B+E, A, 0, B} and I = 2. Networks a) and c)-e) 
are weakly reversible, and all except f) have deficiency 5 = 0. 



{Ai, A 2 , A 3 , . . .}, and these evolve according to c = f (c). 
If we assume mass action kinetics (as we do in most of 
this article), then f(c) = J2 y , y > k y^v' cV (y' ~y)> where y 
and y' denote reactant and product complexes, k y ^ y i is 
the rate constant for the reaction y — > y 1 (and is if the 
reaction does not occur) , and y is a vector of length ns 
giving the stoichiometric coefficient y^ of each species in 
complex y [e.g. if y = 2A!+3A 2 , then y = (2, 3, 0, 0, . . .)]. 
The expression c y is shorthand for Ilj=i ^a - 3 • 

This f(c) can be rewritten in a form that will prove 
useful once we add stochastic effects: 



f (c) = YA(fc)*(c) 



(1) 



Here, Y is an rig x nc matrix whose columns are the 
stoichiometric vectors y of the tlq complexes. A is an 
Tic x riQ matrix with entries (indexed by the names of 
the complexes) 



ky — *y f i V 7^ U 

J2 y " ky^y"> y = y 



(2) 



The argument k indicates that A(fc) depends on the 
rate constants vy i as well as on the network's topol- 
ogy. Finally, ^(c) is a vector of length uq with elements 

%(c) = c y. 

The DRE c = f(c) gives a good account of the evo- 
lution of mean concentrations in macroscopic systems, 
but it hides the fact that individual reaction events oc- 
cur stochastically. A kinetic equation that incorporates 
this effect becomes important when fluctuations are large 
and is also essential to any statistical-mechanical treat- 
ment of reaction networks. This equation is the chemical 
master equation (CME) [6j; it specifies how the proba- 
bility -P(m, t) that there are m = (rriA 1 , tja 2 i m A 3 , ■ ■ ■) 
molecules of the various species at time t changes in time. 



The CME takes the form 



dP(m,t) 
dt 



— K y^y' 
y,y' 



(m + y-yQ! , 

— — — P m + y - y 

(m-y' ! 



m! 



(m - y)! 
>VP(m) , 



P(m) 



(3) 



where, for a vector v, v! = f^Wj!, with the product taken 
over all components Vj, and the factors m!/(m — y)! are 
the generalization to discrete molecule numbers of c y in 
the DRE. The rate constants K y ^ y > differ from the con- 
stants k y ^ y i in the DRE for the same network by a factor 

of V~ 1+ ^-- j VA i , where V is the system volume. 

Because Eq. [3] is a master equation, it must have a 
SSD Pss( m ) satisfying WPss( m ) = 0; with trivial ex- 
ceptions, Pss is unique Q. For a dilute chemical system 
in equilibrium with a heat bath, we know that Pss must 
take the Poisson form 



Pss(m) 



i c m 

Z m! 



= P(m;C) • 



(4) 



The individual components (a of C include factors of the 
volume V and of the thermal wavelength and may incor- 
porate a partition sum over internal molecular degrees of 
freedom. In general, a reaction network has conservation 
laws of the form y • m = constant, and the constant Z is 
chosen so that Pss is normalized for given values of the 
conserved quantities. 

We now ask whether any far-from-equilibrium reaction 
networks have a SSD of the Poisson form Q. We can 
answer this question simply by substituting a Poisson 
distribution into the equation WPss( m ) = 0. One finds 
that 



WP(m;C) cx 

y,y' 



*m+y-y 



(m-y')! (m-y)! 
! 



(5) 



m! 



[m-y 



A(k)*(C) , 



where (, ) is the standard inner product in R nc and u y 
is a unit vector with the component corresponding to y 
equal to 1 and all others 0; the set of w y for all complexes 
y is a basis for K™ c . The second line of Eq. O makes clear 
the connection between the existence of a PSSD and the 
DRE in the form (TJ. 

For a given network topology and choice of the K y ^ y > , 
there is a PSSD if WP(m; Q = for some £ and all m. 
Clearly, it is sufficient that the following condition hold: 



There exists £ such that A («)*(£) = 0. 



(6) 



Because the left side of the inner product in Eq. [5] de- 
pends on m while the right side does not, one might 
guess that it is also necessary that the right side van- 
ish to have WP = for all m. This is indeed the case. 
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In brief, one constructs a set of values of m such that 
the corresponding vectors on the left side of the inner 
product in (JSJ) span M. nc : beginning with an m such that 
m!/(m — y)! vanishes for all but one complex y, one in- 
creases the components of m until the ratio is nonzero 
for two complexes, and so on until the space is spanned. 
(When the network has conserved quantities, the argu- 
ment requires the introduction of a generalized chemical 
potential.) 

The preceding paragraph gives our main result: For a 
reaction network to have a PSSD, it is necessary and suf- 
ficient that ([6|) hold. The same equation appears in the 
theory of DRE's. From Eq. [TJ it is clear that if there is 
a concentration vector c* such that A(fc)5'(c*) = 0, this 
c* is a fixed point of the dynamical system c = f(c). 
It turns out that one can say considerably more p}- 
Specifically, if there exists such a c*, it is the unique 
fixed point (for given values of any conserved quanti- 
ties), and it comes equipped with a Lyapunov function 
£(c) - ££i ca 3 In(c Aj /<% ) - (c Aj - c A .). This strongly 
resembles the free energy of a mixture of ideal gases 
with concentrations c A • Thus, strikingly, those reac- 
tion networks that have an equilibrium-like Poisson SSD 
when stochastic effects are taken into account are pre- 
cisely those that minimize a free-energy-like function in 
the DRE limit. 

The theory of Feinberg and coworkers Q also gives 
an easy way to identify many networks that satisfy (|6j). 
Every network can be assigned a topological index, the 
deficiency 8 = nc — I — s. Here, s is the dimension of the 
stoichiometric subspace S of R™ s spanned by all vectors 
y — y' with K«— ^ 0, and I is the number of linkage 
classes. Two complexes y and y' are in the same linkage 
class if there is a chain of reactions y — > y% — ► y-x — ► . . . — » 
y' linking either y to y' or y' to y. A network is said to be 
weakly reversible if whenever there is such a chain from 
y to y' there is also one from y' to y (Fig. [I} . 

References 0] show that ([6]) holds for any weakly re- 
versible network with (5 = 0; such networks thus all have 
a PSSD. What sort of networks are these? Fig.QJi-e gives 
some examples. Since s < ns is a measure of the number 
of independent species in the network, we can interpret 
the requirement 8 = as a sort of sparseness condition: 
In a deficiency zero network, the average molecule cannot 
react with too many different sets of partners — that is, be 
in too many different complexes. It is this sparseness that 
gives these networks their equilibrium-like character. For 
an equilibrium network, we can define the standard free 
energy difference between two species by summing the 
log ratio of forward to backward rates along a reaction 
path that converts one species into the other. One can 
think of rn(£ A . ) m the PSSD for a deficiency zero network 
as being constructed in a roughly similar manner. The 
difference is that in the one case detailed balance guar- 
antees that a consistent £ can be found no matter which 
path we take between species. In the other it is instead 
the sparseness enforced by 8 = that ensures that no 
inconsistencies arise. 



We next turn to some extensions of our central result 
describing the set of networks with a PSSD. The abil- 
ity to satisfy a fluctuation-dissipation theorem (FDT) 
is often taken as a hallmark of equilibrium-like behav- 
ior [H, 0). It is less well-appreciated that any system 
governed by a master equation obeys a similar linear re- 
sponse relation In general, if the SSD Pss( m !A i ) 
depends on a parameter the observable conjugate 
to fx is M (m) = (fe B T)91n[Pss(m;/z)]/<9/z. That is, 
for any other observable O'(m), the correlation func- 
tion {O'(t)Ofj,(0)) determines the linear response of O' 
to a time-dependent n(t), just as in the FDT. A driven 
system can thus be said to obey the FDT if the con- 
jugate parameter-observable pairs (^,O m ) are those ex- 
pected from equilibrium reasoning. For a system with a 
PSSD, this is easily shown. We focus here on the simple 
case where fi plays the role of a chemical potential, but 
the generalization to other parameters is straightforward. 

To model contact with a particle reservoir, we add to 
any network the reactions 4> ^ y^ with rates k+ to create 
and K- to destroy particles, where y^ is a complex with 
only 1 species. A large enough reservoir will remain in 
equilibrium even when brought into contact with a driven 
reaction network; its chemical potential then turns out 
to be \i = ln(K + / K- ) plus constants that depend only on 
the system volume and y^s composition. One can readily 
show that if the original network with stoichiometric sub- 
space S had a PSSD, the new system in contact with the 
reservoir also has a PSSD provided y M ^ S. The param- 
eters C,^ characterizing the new PSSD are related to the 
parameters C of the old PSSD by C M = Ce u - Here the j th 
component of £e u is Qe Uj , and u oc (/i + constants), 
with y^ 1 - the projection of y^ onto the orthogonal com- 
plement S 1 " of S. Substituting £^ into the expression (|U) 
defining a PSSD, we see that the observable conjugate to 
the chemical potential [i is y^ 1 - • m, just as in an equilib- 
rium system. The driven system thus appears to satisfy 
the FDT. 

A related phenomenon occurs if two reaction ves- 
sels are brought into contact and allowed to exchange 
molecules of species A (Fig. [2]). Define the (generalized) 
chemical potential fi of the molecules in a given vessel by 
adding the reactions <f> A as in the previous paragraph 
and choosing e M oc k + / k_ so that the average molecule 
number of each species remains unchanged. One can 
show that the A molecules in the two vessels have the 
same chemical potential at steady state, again in perfect 
analogy to equilibrium thermodynamics. 

So far, we have dealt only with reactions in well-stirred 
systems. It turns out that our results remain true even 
without stirring. Specifically, consider a system in which 
the same reaction network is replicated at each site on 
a lattice, and a molecule A^ can hop to any neighbor- 
ing site with a rate Xj . If the reaction network at each 
site in isolation has a PSSD, the same is true of the cou- 
pled system on the lattice, whose SSD takes the form 
Pss(m 1 ,m 2 ,m 3 ,...) = (1/^)C E ° m ° AIL m " ! ), wher e 
m Q gives the number of each species on lattice site a. To 
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FIG. 2: Two reaction networks separated by a semi- 
permeable membrane through which they exchange A 
molecules. At steady state, A has the same (generalized) 
chemical potential on both sides. 

verify this result, simply observe that the terms propor- 
tional to each Xj in WPss separately sum to zero. 

Finally, we note that, as in the deterministic case [l^ ]. 
our results extend to certain models with non-mass- 
action kinetics: If the contribution toa j !/( to a j — UAjV- 
of the species Aj to the reaction rate of the complex 
y is everywhere replaced by rii=m A .-y A .+i w &j (i), for 
some function WAj , then the PSSD generalizes to a prod- 
uct over species of the form FJ A 7V\. (mx ) ■ Such mod- 
els arise as approximate descriptions of networks with 
enzymatically catalyzed reactions; the product form of 
the SSD implies that noise in one species does not feed 
through the network to other species downstream [l3|. 
Our results thus extend [l3| to reactions involving more 
than one reactant. 

In sum, we have given the first characterization of the 
set of chemical reaction networks whose stationary state 
is a Poisson distribution, as it would be for a system in 
equilibrium. It turns out that a large class of networks 
have this property, including all those that are weakly 
reversible and are sufficiently sparse, in the sense that 
they have 6 = 0. Whether 6 vanishes depends only on 
which reactant species can be turned into which prod- 
uct species; a network with the appropriate topology will 
have a PSSD for any values of its rate constants. This 
contrasts with the equilibrium case, where varying a rate 
generally drives a network out of equilibrium without any 
change in its topology. 

The significance of these results is severalfold. In light 



of the tremendous interest in noise in biochemical sys- 
tems, there is a need to strengthen the theoretical foun- 
dations of the subject. In particular, there are very few 
exact results for non-trivial networks [TjJ ; in most cases, 
the only techniques available are simulations or expan- 
sions about the DRE. This paper substantially increases 
our arsenal of solvable models. It appears that few func- 
tional biological networks have 6 = 0; indeed, it is possi- 
ble that biology avoids such architectures in part because 
their PSSD enforces a variance in molecule number of or- 
der the mean. Nonetheless, one can imagine many uses 
for our results. The most obvious is as a starting point for 
a perturbative treatment of networks near those with a 
PSSD; this would represent a set of approximations inde- 
pendent of expansions in noise strength, and thus would 
complement such standard approaches. 

As important as any biological applications is the new 
perspective our work gives on the physics of driven sys- 
tems. There has long been interest in equilibrium-like 
descriptions of non-equilibrium phenomena, and chemi- 
cal networks with a PSSD represent a novel set of exam- 
ples along these lines. In more extensively-studied gran- 
ular systems, a common intuition is that a thermal de- 
scription should apply when one energy scale dominates 
and can be used to define an effective temperature. One 
can think of this energy scale as arising from a balance 
between energy input and dissipation. Chemical reac- 
tion networks suggest a generalization of this scenario. 
A thermal analogy holds when the number of interac- 
tions available to a given species, and thus the number 
of channels to inject or dissipate energy, is limited. This 
defines an energy scale for each species; differences be- 
tween these scales are accommodated by our freedom to 
choose ln(CA,0) which plays the role of an internal free 
energy per molecule, for each A.,-. It is tempting to spec- 
ulate that the equilibrium analogy could be extended to 
other driven systems with a similar flexibility to define 
different energy scales for different components. 

I am grateful to Eduardo Sontag for helpful discus- 
sions. This work was funded in part by NSF grant 
PHY05-51164 to the KITP and by FOM/NWO. 



[1] R.P. Ojha, P.-A. Lemieux, P.K. Dixon, A.J. Liu, and 
D.J. Durian, Nature 427, 521 (2004); A.R. Abate and 
D.J. Durian, Phys. Rev. E 72, 031305 (2005). 

[2] G.W. Baxter and J.S. Olafsen, Nature 425, 680 (2003). 

[3] C. Song, P. Wang, and HA. Makse, Proc. Nat'l. Acad. 
Set. USA 102, 2299 (2005). 

[4] I.K. Ono et al, Phys. Rev. Lett. 89, 095703 (2002); 
C.S. O'Hern, A.J. Liu, and S.R. Nagel, ibid. 93, 165702 
(2004); N. Xu and C.S. O'Hern, ibid. 94, 055701 (2005). 

[5] Y. Shokef, G. Bunin, and D. Levine, Phys. Rev. E 73, 
046132 (2006). 

[6] M. Kaern, T.C. Elston, W.J. Blake, and J.J. Collins, 
Nature Rev. Genetics 6, 451 (2005); M.S. Samoilov and 
A.P. Arkin, Nature Biotech. 24, 1235 (2006); P.S. Swain 



and A.Longtin, Chaos 16, 026101 (2006). 

[7] M. Feinberg, Chem. Eng. Set. 42, 2229 (1987); M. Fein- 
berg, Arch. Rat. Mech. Anal. 132, 311 (1995). 

[8] For a pedagogical treatment of the material 
in |7j], see M. Feinberg, unpublished lecture notes, 
ht t p : / / www .che.eng.ohio-state.edu/~feinberg 
/LecturesOnReactionNetworks/; J. Gu- 

nawardena, unpublished lecture notes, 

|http://www.jeremy-gunawardena.com /papers /crnt.pdf[ 

[9] N. van Kampen, Stochastic Processes in Physics and 
Chemistry (2nd ed.). North Holland, 2001. 
[10] R. Graham, Z. Phys. B 26, 397 (1977). 
[11] J.E. Hornos et a/., Phys. Rev. E 72, 051907 (2005); T. 
Fournier et al., Bioinformatics 23, 3185 (2007) 



■5 



[12] E.D. Sontag, IEEE Trans. Auto. Contr. 46, 1028 (2001). 9224 (2007). 

[13] E. Levine and T. Hwa, Proc. Nat'l. Acad. Set. USA 104, 



